library(here)
library(plotly)
theme_set(theme_minimal())
source(here("02_scripts/01_fun_data.R"))
t_owsa_ini_pop <- read.csv(file = here("01_data/ows_tbl_ini_pop.csv")) %>% select(-X)
ows_tbl_prob <- read.csv(, file = here("01_data/ows_tbl_prob.csv")) %>% select(-X)# m_owsa_ini_pop <- matrix(0,
# nrow = nrow(ini_states_tbl_org),
# ncol = 8,
# dimnames = list(
# rows = paste0("ini_pop_",v_state_names),
# cols = c("input_low", "input_high",
# "costs_low", "costs_high",
# "deaths_low", "deaths_high",
# "od_deaths_low", "od_deaths_high")
# ))
#
# for (state in v_state_names){
#
# t_ini_pop_temp <- ini_states_tbl_org %>%
# mutate(basevalue = ifelse(`Variable ` == state, lb, basevalue))
#
# mod <- markov_mod(num_cycles = n_cycles,
# vec_state_names = v_state_names,
# vec_m_0 = t_ini_pop_temp$basevalue,
# mat_P = NA,
# array_P = arr_P,
# vec_cost_states = v_cost_states,
# disc_fac = 0.03, t_0_cost = 0, inc = pop_increase_tbl_new$increase)
#
# deaths_lb <- mod$m_M[181, "BO_DEATH"]
# od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# costs_lb <- mod$total_net_present_cost
#
# t_ini_pop_temp <- ini_states_tbl_org %>%
# mutate(basevalue = ifelse(`Variable ` == state, ub, basevalue))
#
# mod <- markov_mod(num_cycles = n_cycles,
# vec_state_names = v_state_names,
# vec_m_0 = t_ini_pop_temp$basevalue,
# mat_P = NA,
# array_P = arr_P,
# vec_cost_states = v_cost_states,
# disc_fac = 0.03, t_0_cost = 0, inc = pop_increase_tbl_new$increase)
#
# deaths_ub <- mod$m_M[181, "BO_DEATH"]
# od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# costs_ub <- mod$total_net_present_cost
#
# m_owsa_ini_pop[paste0("ini_pop_",state), ] <- c(ini_states_tbl_org$lb[ini_states_tbl_org$`Variable ` == state],
# ini_states_tbl_org$ub[ini_states_tbl_org$`Variable ` == state],
# costs_lb, costs_ub,
# deaths_lb, deaths_ub,
# od_deaths_lb, od_deaths_ub)
# }
#
#
# t_owsa_ini_pop <- as_tibble(m_owsa_ini_pop) %>%
# mutate(costs_range = abs(costs_high - costs_low),
# deaths_range = abs(deaths_high - deaths_low),
# od_deaths_range = abs(od_deaths_high - od_deaths_low),
# name = rownames(m_owsa_ini_pop))t_owsa_ini_pop_costs <- t_owsa_ini_pop %>%
select(costs_low, costs_high, name)
t_owsa_ini_pop_costs$name <- reorder(t_owsa_ini_pop_costs$name, t_owsa_ini_pop$costs_range)
t_owsa_ini_pop_costs <- t_owsa_ini_pop_costs %>%
pivot_longer(1:2, names_to = "group", values_to = "costs") %>%
mutate(costs_new = (costs - mod_basecase$total_net_present_cost)) %>%
filter(costs_new != 0) %>%
mutate(group = ifelse(group == "costs_high", "+25%",
ifelse(group == "costs_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_costs %>% filter(abs(costs_new) >= 500000000),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))t_owsa_ini_pop_deaths <- t_owsa_ini_pop %>%
select(deaths_low, deaths_high, name)
t_owsa_ini_pop_deaths$name <- reorder(t_owsa_ini_pop_deaths$name, t_owsa_ini_pop$deaths_range)
t_owsa_ini_pop_deaths <- t_owsa_ini_pop_deaths %>%
pivot_longer(1:2, names_to = "group", values_to = "deaths") %>%
mutate(deaths_new = (deaths - mod_basecase$m_M[181, "BO_DEATH"])) %>%
filter(deaths_new != 0) %>%
mutate(group = ifelse(group == "deaths_high", "+25%",
ifelse(group == "deaths_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_deaths %>% filter(abs(deaths_new) >= 5000),
aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Deaths from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use",
"ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop %>%
select(od_deaths_low, od_deaths_high, name)
t_owsa_ini_pop_oddeaths$name <- reorder(t_owsa_ini_pop_oddeaths$name, t_owsa_ini_pop$od_deaths_range)
t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop_oddeaths %>%
pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>%
mutate(oddeaths_new = (oddeaths - (mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths))) %>%
filter(oddeaths_new != 0) %>%
mutate(group = ifelse(group == "od_deaths_high", "+25%",
ifelse(group == "od_deaths_low", "-25%", group)))
ggplotly(ggplot() +
geom_bar(data = t_owsa_ini_pop_oddeaths %>% filter(abs(oddeaths_new) >= 100),
aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in opioid-related overdose deaths from the basecase value") + xlab("") +
labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
"ini_pop_BN_PN" = "Pain free, no use",
"ini_pop_BN_ACUTE" = "Acute pain, no use",
"ini_pop_BN_CANCER" = "Cancer, no use",
"ini_pop_BN_OTHER" = "Other, no use",
"ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
"ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
"ini_pop_BPO_CANCER" = "Cancer, Rx use",
"ini_pop_BPO_OTHER" = "Other, Rx use",
"ini_pop_BPO_MISUSE" = "Rx opioid misuse",
"ini_pop_BI_ILLICIT" = "Illicit opioid use",
"ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use",
"ini_pop_BS_OAT_MAINT" = "OAT maintenance")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))# # gen_arr <- function(m_P_temp) {
# # arr_P_temp <- array(0, dim=c(n_states, n_states, n_cycles),
# # dimnames = list(v_state_names, v_state_names, 1:n_cycles))
# #
# # p_od_illicit_temp <- m_P_temp["BO_OD_ILLICIT", "BI_ILLICIT"] # 1-rest for transitions from OD
# # p_illicit_illicit_temp <- m_P_temp["BI_ILLICIT", "BI_ILLICIT"] # 1-rest for transitions from illicit
# #
# # for (t in 1:n_cycles) {
# #
# # m_P_this_cycle <- m_P_temp
# #
# # if(t %in% year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year %in% c(2015, 2016, 2017)]){
# #
# # arr_P_temp[,,t] <- m_P_this_cycle
# #
# # } else if (t %in% year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year %in% c(2018, 2019, 2020)]) {
# #
# # m_P_this_cycle["BO_OD_ILLICIT","BO_OD_DEATH"] <- m_P_temp["BO_OD_ILLICIT", "BO_OD_DEATH"] * (1.0055)^(t-36)
# # m_P_this_cycle["BO_OD_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BO_OD_ILLICIT"] - p_od_illicit_temp)
# #
# # m_P_this_cycle["BI_ILLICIT", "BO_OD_ILLICIT"] <- m_P_temp["BI_ILLICIT", "BO_OD_ILLICIT"] * (1.0055)^(t-36)
# # m_P_this_cycle["BI_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BI_ILLICIT"] - p_illicit_illicit_temp)
# #
# # } else {
# #
# # m_P_this_cycle["BO_OD_ILLICIT","BO_OD_DEATH"] <- m_P_temp["BO_OD_ILLICIT", "BO_OD_DEATH"] * (1.0055)^(36)
# # m_P_this_cycle["BO_OD_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BO_OD_ILLICIT"] - p_od_illicit_temp)
# #
# # m_P_this_cycle["BI_ILLICIT", "BO_OD_ILLICIT"] <- m_P_temp["BI_ILLICIT", "BO_OD_ILLICIT"] * (1.0055)^(36)
# # m_P_this_cycle["BI_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BI_ILLICIT"] - p_illicit_illicit_temp)
# #
# # }
#
# # arr_P_temp[,,t] <- m_P_this_cycle
# # }
# # return(arr_P_temp)
# # }
# #
# # trans_prob_tbl_new <- trans_prob_tbl_new %>%
# full_join(., trans_prob_tbl %>%
# # select(var_from, var_to, lb, ub), by = c("var_from", "var_to")) %>%
# # replace(is.na(.), 0)# trans_prob_tbl_ows <- cbind.data.frame(rep(v_state_names, each = n_states),
# rep(v_state_names, n_states))
#
# names(trans_prob_tbl_ows) <- c("var_from", "var_to")
# names(trans_prob_tbl) <- c("var_from", "var_to",
# "basevalue", "lb", "ub", "range")
#
# trans_prob_tbl_ows <- trans_prob_tbl_ows %>%
# full_join(., trans_prob_tbl %>%
# select(var_from, var_to, basevalue, lb, ub), by = c("var_from", "var_to")) %>%
# replace(is.na(.), 0)
#
# v_var1 <- unique(trans_prob_tbl_ows$var_from[!(trans_prob_tbl_ows$var_from %in% c("BO_OD_DEATH", "BO_DEATH"))])
# v_var2 <- trans_prob_tbl_ows$var_to[trans_prob_tbl_ows$basevalue == 999]
#
#
# dsa_fun <- function(var1, var2) {
#
# m_temp <- matrix(0,
# nrow = nrow(trans_prob_tbl_ows[trans_prob_tbl_ows$var_from == var1, ]),
# ncol = 8 + 4,
# dimnames = list(
# rows = paste0("p_", trans_prob_tbl_ows$var_to[trans_prob_tbl_ows$var_from == var1]),
# cols = c("input_bc", "input_low", "input_high",
# "costs_low", "costs_high", "cost_bc",
# "deaths_low", "deaths_high", "deaths_bc",
# "od_deaths_low", "od_deaths_high", "od_deaths_bc")
# ))
#
# for(state_to in unique(trans_prob_tbl_ows$var_to)){
#
# t_p_temp_full <- trans_prob_tbl_ows %>%
# mutate(basevalue = ifelse(var_to == state_to &
# var_from == var1 &
# var_to != var2, lb, basevalue))
#
# m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# nrow = n_states, ncol = n_states,
# dimnames = list(v_state_names, v_state_names))
#
# for (i in 1:n_states){
# for (j in 1:n_states){
# m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# }
# }
#
# arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#
# mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#
# deaths_lb <- mod$m_M[181, "BO_DEATH"]
# od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# costs_lb <- mod$total_net_present_cost
#
# t_p_temp_full <- trans_prob_tbl_ows %>%
# mutate(basevalue = ifelse(var_to == state_to &
# var_from == var1 &
# var_to != var2, ub, basevalue))
#
# m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# nrow = n_states, ncol = n_states,
# dimnames = list(v_state_names, v_state_names))
#
#
# for (i in 1:n_states){
# for (j in 1:n_states){
# m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# }
# }
#
# arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#
# mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#
#
# deaths_ub <- mod$m_M[181, "BO_DEATH"]
# od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# costs_ub <- mod$total_net_present_cost
#
#
# t_p_temp_full <- trans_prob_tbl_ows
#
# m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# nrow = n_states, ncol = n_states,
# dimnames = list(v_state_names, v_state_names))
#
#
# for (i in 1:n_states){
# for (j in 1:n_states){
# m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# }
# }
#
# arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#
# mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#
#
# deaths_bc_ows <- mod$m_M[181, "BO_DEATH"]
# od_deaths_bc_ows <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# costs_bc_ows <- mod$total_net_present_cost
#
#
#
#
# m_temp[paste0("p_", state_to), ] <- c(trans_prob_tbl_ows$basevalue[trans_prob_tbl_ows$var_from == var1 &
# trans_prob_tbl_ows$var_to == state_to],
# trans_prob_tbl_ows$lb[trans_prob_tbl_ows$var_from == var1 &
# trans_prob_tbl_ows$var_to == state_to],
# trans_prob_tbl_ows$ub[trans_prob_tbl_ows$var_from == var1 &
# trans_prob_tbl_ows$var_to == state_to],
# costs_lb, costs_ub,costs_bc_ows,
# deaths_lb, deaths_ub,deaths_bc_ows,
# od_deaths_lb, od_deaths_ub, od_deaths_bc_ows)
# }
#
# return(as_tibble(m_temp) %>%
# mutate(group = var1,
# cost_range = abs(costs_high - costs_low),
# deaths_range = abs(deaths_high - deaths_low),
# od_deaths_range = abs(od_deaths_high - od_deaths_low),
# name = rownames(m_temp)))
#
# }# # basecase for ows
#
# # m_P_bc <- matrix(trans_prob_tbl_new$basevalue, byrow = T,
# # nrow = n_states, ncol = n_states,
# # dimnames = list(v_state_names, v_state_names))
# #
# # for (i in 1:n_states){
# # for (j in 1:n_states){
# # m_P_bc[i, j] <- ifelse(m_P_bc[i, j] == 999, 1 - (rowSums(m_P_bc)[i] - 999), m_P_bc[i, j])
# # }
# # }
#
# #
# #
# #
# # modbc <- markov_mod(array_P = NA, mat_P = m_P_bc)
# # deaths_bc <- modbc$m_M[181, "BO_DEATH"]
# # od_deaths_bc <- modbc$m_M[181, "BO_OD_DEATH"] + modbc$extra_od_deaths
# # costs_bc <- modbc$total_net_present_cost
# #
# #
# # dsa_fun <- function(var1, var2) {
# #
# # m_temp <- matrix(0,
# # nrow = nrow(trans_prob_tbl_new[trans_prob_tbl_new$var_from == var1, ]),
# # ncol = 8 + 4,
# # dimnames = list(
# # rows = paste0("p_", trans_prob_tbl_new$var_to[trans_prob_tbl_new$var_from == var1]),
# # cols = c("input_bc", "input_low", "input_high",
# # "costs_low", "costs_high", "cost_bc",
# # "deaths_low", "deaths_high", "deaths_bc",
# # "od_deaths_low", "od_deaths_high", "od_deaths_bc")
# # ))
# #
# # for(state_to in unique(trans_prob_tbl_new$var_to)){
# #
# # t_p_temp_full <- trans_prob_tbl_new %>%
# # mutate(basevalue = ifelse(var_to == state_to &
# # var_from == var1 &
# # var_to != var2, lb, basevalue))
# #
# # m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# # nrow = n_states, ncol = n_states,
# # dimnames = list(v_state_names, v_state_names))
# #
# # for (i in 1:n_states){
# # for (j in 1:n_states){
# # m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# # }
# # }
# #
# # mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #
# # deaths_lb <- mod$m_M[181, "BO_DEATH"]
# # od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# # costs_lb <- mod$total_net_present_cost
# #
# #
# # t_p_temp_full <- trans_prob_tbl_new %>%
# # mutate(basevalue = ifelse(var_to == state_to &
# # var_from == var1 &
# # var_to != var2, ub, basevalue))
# #
# # m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# # nrow = n_states, ncol = n_states,
# # dimnames = list(v_state_names, v_state_names))
# #
# #
# # for (i in 1:n_states){
# # for (j in 1:n_states){
# # m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# # }
# # }
# #
# # mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #
# #
# # deaths_ub <- mod$m_M[181, "BO_DEATH"]
# # od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# # costs_ub <- mod$total_net_present_cost
# #
# #
# # t_p_temp_full <- trans_prob_tbl_new
# #
# # m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# # nrow = n_states, ncol = n_states,
# # dimnames = list(v_state_names, v_state_names))
# #
# # for (i in 1:n_states){
# # for (j in 1:n_states){
# # m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# # }
# # }
# #
# # mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #
# # deaths_bc_tmp <- mod$m_M[181, "BO_DEATH"]
# # od_deaths_bc_tmp <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# # costs_bc_tmp <- mod$total_net_present_cost
# #
# #
# # m_temp[paste0("p_", state_to), ] <- c(trans_prob_tbl_new$basevalue[trans_prob_tbl_new$var_from == var1 &
# # trans_prob_tbl_new$var_to == state_to],
# # trans_prob_tbl_new$lb[trans_prob_tbl_new$var_from == var1 &
# # trans_prob_tbl_new$var_to == state_to],
# # trans_prob_tbl_new$ub[trans_prob_tbl_new$var_from == var1 &
# # trans_prob_tbl_new$var_to == state_to],
# # costs_lb, costs_ub,costs_bc_tmp,
# # deaths_lb, deaths_ub,deaths_bc_tmp,
# # od_deaths_lb, od_deaths_ub, od_deaths_bc_tmp)
# # }
# #
# # return(as_tibble(m_temp) %>%
# # mutate(group = var1,
# # cost_range = abs(costs_high - costs_low),
# # deaths_range = abs(deaths_high - deaths_low),
# # od_deaths_range = abs(od_deaths_high - od_deaths_low),
# # name = rownames(m_temp)))
# #
# # }
# ## ows_tbl_prob <- data.frame()
#
# for (i in 1:length(v_var1)){
# ows_tbl_prob <- rbind.data.frame(ows_tbl_prob, dsa_fun(v_var1[i], v_var2[i]))
# }
#
# write.csv(ows_tbl_prob, file = "../01_data/ows_tbl_prob.csv")ows_tbl_prob_costs <- ows_tbl_prob %>%
select(costs_low, costs_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_costs$name <- reorder(ows_tbl_prob_costs$name, ows_tbl_prob$cost_range)
ows_tbl_prob_costs <- ows_tbl_prob_costs %>%
pivot_longer(1:2, names_to = "group", values_to = "costs") %>%
mutate(costs_new = (costs - as.numeric(mod_basecase$total_net_present_cost)),
costs_new_per = (costs_new/as.numeric(mod_basecase$total_net_present_cost)) * 100) %>%
filter(costs_new != 0) %>%
mutate(group = ifelse(group == "costs_high", "+25%",
ifelse(group == "costs_low", "-25%", group)))
name_cost_1 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 0.5) %>% select(name)
name_cost_2 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 1.5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_1$name),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_2$name),
aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in Costs from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))ows_tbl_prob_deaths <- ows_tbl_prob %>%
select(deaths_low, deaths_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_deaths$name <- reorder(ows_tbl_prob_deaths$name, ows_tbl_prob$deaths_range)
ows_tbl_prob_deaths <- ows_tbl_prob_deaths %>%
pivot_longer(1:2, names_to = "group", values_to = "deaths") %>%
mutate(deaths_new = (deaths - as.numeric(mod_basecase$m_M[181, "BO_DEATH"])),
deaths_new_per = (deaths_new/as.numeric(mod_basecase$m_M[181, "BO_DEATH"])) * 100) %>%
filter(deaths_new != 0) %>%
mutate(group = ifelse(group == "deaths_high", "+25%",
ifelse(group == "deaths_low", "-25%", group)))
name_death_1 <- ows_tbl_prob_deaths %>% filter(abs(deaths_new_per) >= 0.5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_deaths %>% filter(name %in% name_death_1$name),
aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in deaths from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))oddeaths_bc <- as.numeric(mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths)
ows_tbl_prob_oddeaths <- ows_tbl_prob %>%
select(od_deaths_low, od_deaths_high, name, group) %>%
mutate(name = paste(name, group, sep = "_")) %>%
rename(grp = group)
ows_tbl_prob_oddeaths$name <- reorder(ows_tbl_prob_oddeaths$name, ows_tbl_prob$od_deaths_range)
ows_tbl_prob_oddeaths <- ows_tbl_prob_oddeaths %>%
pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>%
mutate(oddeaths_new = (oddeaths - oddeaths_bc),
oddeaths_new_per = (oddeaths_new/oddeaths_bc) * 100) %>%
filter(oddeaths_new != 0) %>%
mutate(group = ifelse(group == "od_deaths_high", "+25%",
ifelse(group == "od_deaths_low", "-25%", group)))
name_oddeath_1 <- ows_tbl_prob_oddeaths %>% filter(abs(oddeaths_new_per) >= 5) %>% select(name)
ggplotly(ggplot() +
geom_bar(data = ows_tbl_prob_oddeaths %>% filter(name %in% name_oddeath_1$name),
aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
ylab("Difference in opioid-related overdose death from the basecase value") + xlab("") +
labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
"p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
"p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
"p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
"p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
"p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
"p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
"p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
"p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
"p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
"p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
"p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
"p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
"p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
"p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
"p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
"p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
"p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
"p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
"p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
"p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
"p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",
"p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
"p_BO_DEATH_BN_OTHER" = "Other no use to death",
"p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
"p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death",
"p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
"p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
"p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
"p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay",
"p_BO_OD_DEATH_BO_OD_ILLICIT" = "Illicit opioid overdose to death",
"p_BO_OD_ILLICIT_BI_ILLICIT" = "Illicit opioid use to Illicit opioid overdose",
"p_BO_OD_RX_BPO_MISUSE" = "Rx opioid misuse to Rx opioid overdose",
"p_BR_OD_ILLICIT_BR_ILLICIT" = "R: illicit opioid use to R: illicit opioid overdose",
"p_BO_OD_DEATH_BR_OD_ILLICIT" = "R: Illicit opioid overdose to death",
"p_BO_OD_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to Illicit opioid overdose",
"p_BI_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to illicit opioid use")) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2"))